Sodium LGS

Simulate a sodium LGS beacon

Import the LGSSystem class from the pylgs module:

from pylgs.all import *
import numpy as np
import pymor

Define an LGS system on the D\(_2\) line with pump and repump beams, with all parameters fixed except for the light intensities and the magnetic-field strength:

lgs = LGSSystem(
    'NaD2_Repump', 
    {
        'EllipticityDegrees1': 45.,
        'PolarizationAngleDegrees1': 0,
        'DetuningHz1': 1.0832e9,
        'LaserWidthHz1': 10.0e6,
        'EllipticityDegrees2': 45.,
        'PolarizationAngleDegrees2': 0,
        'DetuningHz2': -6.268e8,
        'LaserWidthHz2': 10.0e6,
        'MagneticZenithDegrees': 45.,
        'MagneticAzimuthDegrees': 45.,
        'SDampingCollisionRatePerS': 4081.63,
        'BeamTransitRatePerS': 131.944,
        'VccRatePerS': 28571.,
        'TemperatureK': 185.,
        'RecoilParameter': 1
    }
)

The model contains equations for 374 density-matrix elements:

lgs.A_ind.range
{Density matrix (range)(374)}

Define sample values for the variable parameters:

params = {'IntensitySI1': 5., 'IntensitySI2': 46., 'BFieldG': 0.5}

Steady state of a cw LGS

Evenly spaced velocity groups

Build a model for the steady state of a cw LGS with 200 evenly spaced velocity groups:

model = lgs.stationary_model(vg=200)

Find the steady-state return flux:

model.total_flux(params).item()
13625.270307926916

Find the density matrix for each of the 200 velocity groups:

sol = model.solve(params)
sol
{Atomic velocity(200) ⨉ Density matrix (source)(374), }

The result is expressed as an XarrayVectorArray object representing 200 vectors each of size 374. The XarrayVectorArray is a wrapper for an xarray DataArray, which can be accessed via the XarrayVectorArray.array attribute:

sol.array
<xarray.DataArray (Atomic velocity: 200, Density matrix (source): 374)> Size: 598kB
array([[ 4.39820989e-07,  1.36790153e-07, -1.36692297e-07, ...,
         4.11349076e-11,  2.15032837e-13,  3.90993999e-14],
       [ 5.25623820e-07,  1.63475922e-07, -1.63358570e-07, ...,
         4.96457650e-11,  2.62098981e-13,  4.73161225e-14],
       [ 6.27016980e-07,  1.95009270e-07, -1.94868909e-07, ...,
         5.98154463e-11,  3.18954046e-13,  5.71688740e-14],
       ...,
       [ 6.25164298e-07,  1.95204690e-07, -1.95037857e-07, ...,
        -6.26906845e-11,  3.50162745e-13,  6.07662656e-14],
       [ 5.24087271e-07,  1.63638177e-07, -1.63498877e-07, ...,
        -5.20072904e-11,  2.87471533e-13,  5.02408015e-14],
       [ 4.38562256e-07,  1.36929700e-07, -1.36813585e-07, ...,
        -4.30715543e-11,  2.35630633e-13,  4.14736854e-14]])
Coordinates:
  * Atomic velocity          (Atomic velocity) float64 2kB -2.985 ... 2.985
  * Density matrix (source)  (Density matrix (source)) <U68 102kB 'ρ<sub>Re, ...

Use the density-matrix solution to find the velocity distribution of flux as an XarrayVectorArray:

flux = model.flux_distribution(sol)
flux
{Atomic velocity(200) ⨉ Transition(1), }

XarrayVectorArray objects have built-in visualization routines. Plot the velocity distribution of return flux:

model.flux_distribution(sol).visualize(xaxis_range=(-.5, .5))

Plot the level population distribution:

model.level_population_distribution(sol).visualize()

Adaptive refinement of velocity groups

Build a steady-state model with each velocity group contributing no more than 2% of the return flux:

model = lgs.adaptive_stationary_model(params, max_weight=0.02)

Find the total flux for this model:

model.total_flux(params).item()
13521.856397884763

Find the density matrix. In this case the adaptive refinement has produced 82 velocity groups:

sol = model.solve(params)
sol
{Atomic velocity(82) ⨉ Density matrix (source)(374), }

The flux distribution shows that this is a more detailed description of the system than the model with a greater number of evenly spaced velocity groups:

model.flux_distribution(sol).visualize(xaxis_range=(-.5, .5))

Dependence on system parameters

Find the steady-state density matrix for 10 different values of repump intensity, with the total pump + repump intensity held constant:

sol = model.solve([
    {'BFieldG': .5},
    {
        'IntensitySI1': np.linspace(0, 25, 10), 
        'IntensitySI2': 50 - np.linspace(0, 25, 10)
    }
])
sol
{Atomic velocity(82) ⨉ Density matrix (source)(374), IntensitySI1(10)}

Plot the total flux to show the optimal value of repump intensity for these conditions:

model.total_flux(sol).visualize()

Create an animation of the flux distribution as the repump intensity is varied:

model.flux_distribution(sol).visualize(xaxis_range=(-.5, .5))

Animate the level population distribution as a function of repump intensity:

model.level_population_distribution(sol).visualize()

Transient dynamics

Transient dynamics from a laser pulse

Build a dynamical model (InstationaryModel) using the velocity groups from the adaptively refined steady-state model:

model = lgs.instationary_model(
    vg=model.data['velocity_groups'], 
    T=1e-5, 
    num_values=100
)

Set tolerances for the BDF solver:

pymor.basic.set_defaults({
    'pylgs.pymor.timestepping.cvode_solver_options.cvode_bdf_atol': 1e-6,
    'pylgs.pymor.timestepping.cvode_solver_options.cvode_bdf_rtol': 1e-4,
})

Solve for the density-matrix evolution. The solution contains the density matrix for each velocity group at 100 time values:

sol = model.solve(params)
sol
{Atomic velocity(82) ⨉ Density matrix (source)(374), Time(100)}

Plot the total flux as a function of time:

model.total_flux(sol).visualize()

Animate the flux distribution as a function of time:

model.flux_distribution(sol).visualize()

Animate the level population distribution:

model.level_population_distribution(sol).visualize(xaxis_range=(-1, 1))

Transient dynamics with modulated laser intensity

Solve the InstationaryModel with sinusoidally modulated light intensity:

sol = model.solve({
    'IntensitySI1': "5.*sin(1e6*t)**2", 
    'IntensitySI2': "46.*sin(1e6*t)**2",
    'BFieldG': .5
})

The total flux as a function of time:

model.total_flux(sol).visualize()

The animated flux distribution:

model.flux_distribution(sol).visualize(xaxis_range=(-.25, .25))

The animated level population distribution:

model.level_population_distribution(sol).visualize(xaxis_range=(-.25, .25))

Periodic state of a modulated system

The StationaryFloquetModel solves for the periodic state of a modulated system after the transient dynamics have died out.

Build a stationary Floquet model for six Fourier harmonics of a system with sinusoidally modulated intensity:

model = lgs.stationary_floquet_model(
    {
        'IntensitySI1': "5.*sin(nu*t)**2", 
        'IntensitySI2': "46.*sin(nu*t)**2"
    },
    T=np.linspace(0, 1e-5, 100),
    n_vector=6,
    vg=model.data['velocity_groups']
)

Define values for the modulation frequency and magnetic field:

mu = {'nu': 1e6, 'BFieldG': .5}

Solve for the periodic state. The solution contains 13 Fourier coefficients for each density-matrix element in each velocity group:

sol = model.solve(mu)
sol
{Fourier coefficient(13) ⨉ Atomic velocity(82) ⨉ Density matrix (source)(374), }

The real and imaginary parts of the Fourier coefficients of the time-dependent return flux:

model.total_flux(sol).visualize()

The total return flux reconstructed as a function of time:

model.reconstructed_total_flux(sol, mu).real.visualize()

The level population distribution as a function of time:

model.level_population_distribution(model.reconstruct(sol, mu)).real.visualize()